;+
; NAME:
;       GENERATE_EMAF_LIST
;
; PURPOSE:
;       Generate a list of EMAF sources in the catalog, and create an
;       IDL structure containing relevant information for those
;       sources.
;
; CATEGORY:
;       EMAF
;
; CALLING SEQUENCE:
;       GENERATE_EMAF_LIST
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE  [IDL Save file: ./emaf_paper/bgps_emaf.sav]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  10/23/12, TPEB -- Initial version, coppied and
;                                   significantly culled from
;                                   create_bgps_rb3_struct.pro.
;       Modified: 10/30/12, TPEB -- Added grs_angsep elements to the
;                                   EMAF structure.
;       Modified: 11/01/12, TPEB -- Check for the existance of
;                                   morph_distances.sav and remove --
;                                   that way, any update of the source
;                                   list will cause irdc_morph to be
;                                   re-run.
;       Modified: 11/02/12, TPEB -- Restrict sources to C >= 0.05, and
;                                   I_MIR <= 100 MJy/sr.
;       Modified: 11/29/12, TPEB -- Added SFA structure element to
;                                   encode star-formation activity
;                                   level from Dunham et al (2011,
;                                   ApJ, 731, 90).
;       Modified: 03/11/13, TPEB -- Changed GOTO control statements to
;                                   CONTINUE.
;
;-

PRO GENERATE_EMAF_LIST
  
  ;;===============================================
  ;; Read in information relevant to this routine
  
  ;; Read in the BGPS distance database & velocity structure
  s = read_bgps_csv('bgps_distance_database.csv',csv,/ver)
  restore,'local/bgps_velocity_struct.sav',/ver
  ppbeam = 23.8027594254d       ; # of BGPS pixels per beam
  
  ;; Read in Cara's temperature maps & get useful information
  TMAP_DIR = 'irdc_dist_model/battersby/'
  tmap = readfits(TMAP_DIR+'gaussfit_iter_beta175_temp_l30_conv25.fits',$
                  thdr,/silent)
  tlab = readfits(TMAP_DIR+'gaussfit_iter_beta175_column_l30_'+$
                  'conv25_label.fits',lhdr,/silent)
  terr = readfits(TMAP_DIR+'gaussfit_iter_beta175_temperr_l30_conv25.fits',$
                  ehdr,/silent)
  plot_wcs_axes,thdr,10,lonarr=larr
  lrange = minmax(larr)
  extast,thdr,tastr
  tsz = size(tmap, /DIM)
  
  
  ;;===============================================
  ;; Determine sources w/ dense-gas vlsr & overlap with GLIMPSE
  
  ;; Figure out how many objects have a dense gas composite spectrum
  dgind = WHERE( total(velocity.densegas,1) NE 0, ndg)
  velocity = velocity[dgind]
  
  ;; Check the UBC_PROC files for a list of filenames
  UBC_DIR = './local/ubc_glimpse_proc/'
  ufn = FILE_SEARCH(UBC_DIR+'bgps????.fits',COUNT=n_ubc)
  ubc_cnum = long(strmid(ufn,STRPOS(ufn[0],'bgps')+4,4))
  
  ;; Get the overlap between dense gas tracers and the UBC GLIMPSE proc
  ind = WHERE_ARRAY(velocity.cnum, ubc_cnum, n)
  ufn = ufn[ind]
  
  ;; Work backwards to get UBC with dense gas...
  ind2 = WHERE_ARRAY(ubc_cnum, velocity.cnum, n2)
  velocity = velocity[ind2]
  IF n NE n2 THEN $
     message,"Take a shot of gin, then figure out the size mismatch."
  
  print,'Total sample is '+string(n,format="(I0)")+' objects.'
  
  
  ;;===============================================
  ;; Create a structure to hold the information relevant to selected EMAFs
  
  emaf = replicate( $
         create_struct('name','','cnum',0L,'l',0.,'b',0.,'vlsr',0.,'v_flag',0b,$
                       's40',0.,'s40_err',0.,$
                       'I_mir',0.,'I_min',0.,'C',-1000.,'sc',-1000.,$
                       'sig_imir',0.,'sig_imin',0.,'sig_c',0.,$
                       'higal_td',0.,'higal_te',0,$
                       'grs',0b,'grs_match',0b,'grs_clemdist',0.,'grs_dist',0.,$
                       'grs_angsep',0.,'grs_l',0.,'grs_b',0.,'grs_v',0.,$
                       'grs_21cm','',$
                       'dml',[0.,0.,0.],'dbar',[0.,0.],'p_ratio',0.,'sfa',0,$
                       'n200',0),n)
  
  
  
  
  ;;====================================================================
  ;; Now, loop through each potential source...
  FOR i=0L, n-1 DO BEGIN
     
     ;; Appraise the user of status
     IF ((i+1) MOD 100) EQ 0 THEN message,'Doing object #'+$
                                          string(i+1,format="(I4)")+' of '+$
                                          string(n,format="(I4)"),/inf
     
     ;; Since velocity and ubc are now the same size, loop through
     ;; velocity, since its cnum are in proper order
     j = where(ubc_cnum EQ velocity[i].cnum, nj)
     k = where(s.cnum EQ velocity[i].cnum, nk)
     IF nj EQ 0 THEN message,$
        "Take a shot of rum, then figure out why this object isn't in ubc."
     
     ;; Start filling in the emaf structure from the Bolocat
     emaf[i].name = s[k].name
     emaf[i].cnum = s[k].cnum
     emaf[i].l    = s[k].glon_peak
     emaf[i].b    = s[k].glat_peak
     
     ;; INSERT PROPER FLUX MEASUREMENT HERE -- structure expects Jy
     ;; Bolocat fluxes x 1.5+-0.15
     emaf[i].s40 = s[k].flux_40 * 1.5d
     sig_ss      = (s[k].err_flux_40 > 0.015) * 1.5d
     emaf[i].s40_err  = sqrt(sig_ss*sig_ss + emaf[i].s40*emaf[i].s40*0.0225d)
     
     ;; Insert VLSR from velocity.densegas
     yfit = GAUSSFIT(velocity[i].v_std,velocity[i].densegas, A, NTERMS=3)
     emaf[i].vlsr = A[1]
     emaf[i].v_flag = velocity[i].has_rv_types
     
     
     ;;=============================================================
     ;; Cull based on latitude + location within the L-V diagram:
     ;; Check latitude -- needs to be within safe boundaries of GLIMPSE
     IF (abs(emaf[i].b) GT 1. || emaf[i].l GT 65.25) THEN CONTINUE
     
     ;; Check on the longitude range -- can't use objects w/in 21
     ;;    deg of GC UNLESS part of Molecular Ring structure
     IF (emaf[i].l LT 21.) THEN $
        IF (emaf[i].vlsr GT (10.d/3.d*emaf[i].l + 15.d)) OR $
        (emaf[i].vlsr LT (20.d/9.d*emaf[i].l - 50.d/3.d)) THEN CONTINUE
     
     
     
     ;;=============================================================
     ;; Next, load UBC postage stamp images (created by
     ;; process_glimpse_irdc.pro) and insert mid-infrared data into
     ;; the emaf structure.
     
     scnum = string(s[k].cnum,format="(I04)")
     imir_img = readfits(UBC_DIR+'Imir'+scnum+'.fits',$
                         imir_hdr,/silent)
     nons_img = readfits(UBC_DIR+'nonsmooth'+scnum+'.fits',$
                         nons_hdr,/silent)
     smoo_img = readfits(UBC_DIR+'smooth33arc'+scnum+'.fits',$
                         smoo_hdr,/silent)
     
     ;; For postage-stamp images going off the edge of the GLIMPSE
     ;; map, replace 0's with NANs -- which can more easily be
     ;; dealt with later.
     inan = WHERE(nons_img EQ 0, nnan)
     IF nnan THEN BEGIN
        imir_img[inan] = !values.f_nan
        nons_img[inan] = !values.f_nan
        smoo_img[inan] = !values.f_nan
     ENDIF
     
     ;; TO CALCULATE I_MIR & I_min (with their uncertainties), we use
     ;;    the 40" Bolocat aperture (to align with BGPS flux
     ;;    measurement).
     ;; Extract pixel information using the 40" aperture method
     ;;    from Bolocat
     wtmask  = aperture_wtmask(smoo_hdr, 40., emaf[i].l, emaf[i].b)
     ind     = where(wtmask GT 0.5)
     wtmask2 = aperture_wtmask(imir_hdr, 240., emaf[i].l, emaf[i].b)
     ind2    = where(wtmask2 GT 0)
     
     ;; Calculate fraction of nonsmooth with I >= 200 MJy/sr
     ;;   w/in 2' of peak location
     wtmask3 = aperture_wtmask(imir_hdr, 240., emaf[i].l, emaf[i].b)
     if200 = where((nons_img*wtmask3) GE 200., n200)
     emaf[i].n200 = n200
     
     ;; Calculate the extracted values
     emaf[i].I_min = min( smoo_img[ind], /NAN )
     emaf[i].I_MIR = total(wtmask2*imir_img,/NAN) / total(wtmask2)
     
     ;; Calculate the uncertainties
     emaf[i].sig_imin = stddev( smoo_img[ind], /NAN )
     
     V1 = total(wtmask2,/NAN)
     V2 = total(wtmask2*wtmask2,/NAN)
     emaf[i].sig_imir = sqrt(total(wtmask2[ind2]*$
                                   (imir_img[ind2]-emaf[i].I_mir)^2,/NAN) * $
                             (V1 / (V1*V1 - V2)) )
     
     ;; Calculate contrast
     emaf[i].c = ( 1.d - (emaf[i].I_min / emaf[i].I_mir) )
     emaf[i].sig_C  = sqrt( (emaf[i].i_min*emaf[i].i_min/$
                             emaf[i].i_mir/emaf[i].i_mir)*$
                            (emaf[i].sig_imir*emaf[i].sig_imir/$
                             emaf[i].i_mir/emaf[i].i_mir + $
                             emaf[i].sig_imin*emaf[i].sig_imin/$
                             emaf[i].i_min/emaf[i].i_min) )
     
     emaf[i].sc = 1.d - (min( nons_img[ind], /NAN ) / emaf[i].I_mir)
     ;; Check for correspondence with Cara's temperature maps from Hi-GAL
     AD2XY, emaf[i].l, emaf[i].b, tastr, x, y
     IF (x GE 0) && (x LT tsz[0]) && (y GE 0) && (y LT tsz[1]) THEN BEGIN
        
        ;; Next, check to see if location is within studied region
        IF tlab[x,y] THEN BEGIN
           
           emaf[i].higal_td = tmap[x,y]
           emaf[i].higal_te = terr[x,y]
           
           str = emaf[i].name+'  Td = '+$
                 string(tmap[x,y],terr[x,y],FORMAT="(F4.1,' +- ',F4.1)")+$
                 ' K   LABEL = '+$
                 string(tlab[x,y],FORMAT="(I0)")
           print,str
           
        ENDIF
     ENDIF
     
  ENDFOR
  
  ;; Check for existance of morph_distances.sav AND
  ;;   irdc_morph_pvec.sav, and remove. 
  IF FILE_TEST('./emaf_paper/save_files/morph_distances.sav') THEN $
     spawn,'rm ./emaf_paper/save_files/morph_distances.sav'
  IF FILE_TEST('./emaf_paper/save_files/irdc_morph_pvec.sav') THEN $
     spawn,'rm ./emaf_paper/save_files/irdc_morph_pvec.sav'
  
  ;; Save full list for comparisons, as needed
  save,emaf,filename='./emaf_paper/bgps_emaf_uncut.sav',/ver
  
  
  ;; Check for POSITIVE mid-infrared contrast (C >= 0.05)
  ;; Check for sources with suspicious VLSR ( < -30 km/s)
  ;; Check for I_MIR <= 100 MJy/sr
  ;; Check for I_8um >= 200 MJy/sr w/in 2' (n200 < 10 pixels)
  goodind = WHERE( (emaf.c GT 0.05) AND (emaf.vlsr GE -30.) $
                   AND (emaf.I_MIR LE 100.) AND emaf.n200 LT 10, ngood)
  emaf = emaf[goodind]
  
  save,emaf,filename='./emaf_paper/bgps_emaf_potential.sav',/ver
  print,'Number of potentially useable objects: ',n_elements(emaf)
  
  ;; Check against BY-EYE REJECTION list...
  readcol,'./emaf_paper/rejected_emaf.txt',reject_cnum,format='I'
  good_emaf_ind = MISSING(reject_cnum, emaf.cnum, n_good_emaf)
  print,'# Good: ',n_good_emaf
  emaf = emaf[good_emaf_ind]
  
  ;;====================================================================
  ;; We now have the final list of sources!
  ;; Save the structure
  save,emaf,filename='./emaf_paper/bgps_emaf.sav',/ver
  
  spawn,'cp ./emaf_paper/bgps_emaf.sav /d4/idl_lib/distance-omnibus/emaf/.'
  
  ;; Clear out the memory a liitle
  undefine,velocity,s,ubc,emaf
  
END
